% To do likelihood ratio test
clear all;
clc;
disp('To do likelihood ratio test_MK5');

%load MK2 data (Both tests based on same data)
load PathCountTotal_MK5_TP06YRJfinR5G1000kMR7.mat;
load PathCountAvgProb_MK5_TP06YRJfinR5G1000kMR7.mat;
filename='TP06_Y_RJ_fin_R5_G1000k_MR7_Likelihood Ratio Test_MK4_MK5.xls';

% take the avg data from PathCountAvgProb (MK4).
N5(:,:,:,:,:,:) = PathCountTotal(1,:,:,:,:,:,:);
P5(:,:,:,:,:,:) = PathCountAvgProb(1,:,:,:,:,:,:);
ttlData = sum(N5(:));

% to calculate log (LMK4) hyphothesis
N4(:,:,:,:,:) = sum(N5(1:7,:,:,:,:,:));
P4=zeros(7,7,7,7,8);
LLMK4=0;
for i=1:7
    for j=1:7
        for k=1:7
            for d=1:7
                for t=1:8
                    P4(i,j,k,d,t)=N4(i,j,k,d,t)./sum(N4(i,j,k,d,1:8));
                    if N4(i,j,k,d,t)~=0
                        Y =log(P4(i,j,k,d,t)).* N4(i,j,k,d,t);
                        LLMK4=LLMK4+Y;
                    else
                        LLMK4=LLMK4+0;
                    end
                end
            end
        end
    end
end


% to calculate log (LMK5) hyphothesis
LLMK5=0;
for i=1:7
    for j=1:7
        for k=1:7
            for d=1:7
                for t=1:7
                    for q=1:8
                        if N5(i,j,k,d,t,q)~=0
                            YY =log(P5(i,j,k,d,t,q)).* N5(i,j,k,d,t,q);
                            LLMK5=LLMK5+YY;
                        else
                            LLMK5=LLMK5+0;
                        end
                    end
                end
            end
        end
    end
end


%Likelihood Ratio Test
invalid_Num_N4 = prod(size(find(N4==0)));
invalid_Num_N5 = prod(size(find(N5==0)));
Num_N5=prod(size(P5));
Num_N4=prod(size(P4));
dof = (Num_N5 - invalid_Num_N5) - (Num_N4 - invalid_Num_N4);
%dof = (prod(size(P5))-invalid_Num_N5) - (prod(size(P4))-invalid_Num_N4)% minus the number of zero element in P4 and P5.
%prod(size(P5))-prod(size(P4))-invalid_Num_N4-invalid_Num_N5;
[LRh,LRp,LRstat,cV] = lratiotest(LLMK5,LLMK4,dof);


% PP2(:,:)=P2(2,:,:);
% NN2(:,:)=N2(2,:,:);
% PP3(:,:)=P4(2,2,:,:);
% NN3(:,:)=N4(2,2,:,:);
% sheetlist1=strcat('Prob of ijk_pathR2');
% sheetlist2=strcat('Prob of ijkd_pathR2R2');
% sheetlist3=strcat('ObservationNum N2_pathR2');
% sheetlist4=strcat('ObservationNum N3_pathR2R2');
% sheetlist5=strcat('Parameters');
% xlswrite(filename, PP2, sheetlist1);
% xlswrite(filename, PP3, sheetlist2);
% xlswrite(filename, NN2, sheetlist3);
% xlswrite(filename, NN3, sheetlist4);
% d = {'LRh','LRp','LRstat','cV','dof','LLMK4','LLMK5','ttlData'; LRh, LRp, LRstat,cV,dof,LLMK4,LLMK5,ttlData};
% xlswrite(filename, d, sheetlist5)


